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ABSTRACT 

The kinematic and radiative power of molecular jets is expected to change as a pro- 
tostar undergoes permanent or episodal changes in the rate at which it accretes. We 
study here the consequences of evolving jet power on the spatial and velocity structure, 
as well as the fluxes, of molecular emission from the bipolar outflow. We consider a jet 
of rapidly increasing density and a jet in which the mass input is abruptly cut off. Wc 
perform three dimensional hydrodynamic simulations with atomic and molecular cool- 
ing and chemistry. In this work, highly collimated and sheared jets are assumed. We 
find that position-velocity diagrams, velocity-channel maps and the relative H2 and 
CO fluxes are potentially the best indicators of the evolutionary stage. In particular, 
the velocity width of the CO lines may prove most reliable although the often-quoted 
mass- velocity power-law index is probably not. We demonstrate how the relative H2 1- 
S(l) and CO J=l-0 fluxes evolve and apply this to interpret the phase of several 
outflows. 
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1 INTRODUCTION 

Jets are being increasingly associated with protostars and 
with bipolar molecular outflows. Prominent examples of 
jet-dr i ven collimated outflow s include HH111 ||N^gar_et_alJ 
I l997t [Hartigan et"ai]|200ll). H H211 jGueth fe Guilloteaul 
Il999l: Idiandler fc Richerl l200lh . HH 212 llLee et al.l 1200(3) 
and HH288 JGueth et aljbooili . For many other outflows, 
there is evidence that their momentum is also derived 
from the thrust of jets. For example, transverse shock 
structures beh ind bow shocks ar e interpreted as the sites 
of jet impact (Smi th~et al.ll20o 3). Nevertheless, estimated 
jet thrusts appear insu fficient to drive some outflows 
( Frid lund fc Liseaul Il998t) and ballistic bullets and wide 
winds have been invoked. Alternatively, as investigated 
here, the presently observed jet may not reflect the ancestor 
jet that has previously driven the outflow. 

Our goal is to determine the observable characteristics 
of outflows driven by molecular jets with evolving mass flow 
rates. Jet evolution was considered as an interpretation of 
the systematic decrease in flow velocity with d istance in 
the atomic HH34 outflow iCabrit fc Raeal l200(fl . The de- 
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duced requirements appeared unlikely and an alternative 
precessing-fragmenting scenario was favo ured. A systemat i- 
cally increasing jet speed was simulated bv lLim et alJ {2002). 
Their idea was to explain how molecules can be accelerated 
to high speed without destruction. Here, we choose to evolve 
the mass flow rate rather than the speed since it is clear that 
the mass flow rate must change by many orders of magni- 
tude as a protostar evolves. The taxonomy of protostars is 
defined by the categories Class through Cla ss 3, which 
have been linked to an evolutionary sequence iLadalll987l : 
lAndreet ai1ll993l) . One scenario connects this sequence to 
simultaneous evolutions i n both the accretion rate and the 
outflow rate JSmith 2000 0. Furthermore, certai n protostars 
undergo outbursts iNaravanan fc Walkei)ll99d) . the conse- 
quenc es of which may also show up later in the outflow prop- 
erties jArce fc GoodmarJl200 if) . 

We investigate a limited set of three-dimensional sim- 
ulations that correspond to significant and permanent 
changes to the input mass flux. This complements stud- 
ies of outflow deve lopment generated by non-evolvin ; 
heavy molecular je ts dSuttner et al Jfl"997t iDavis et al 1 ll99; 
IVolker et alJ Il999l) a s well as a range of jet conditions 
lfRo"s^n^^5nu^nl2003l) . We first consider a rapid decrease in 
the mass flux, which we suspected would produce a bullet- 
like outflow. A decrease in mass flux of some variety must 
accompany the transitional stages of a protostar into a star. 
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In addition, we model a jet that increases density after some 
initial propagation, and this model correlates with the oppo- 
site end of a protostar's evolution. During the initial collapse 
of the core into a protostellar disk, a very young jet must ig- 
nite and, therefore, be characterised by an increase in mass 
flux. 

Certain features of the molecular line emission from out- 
flows have been attributed to an evolving jet mass outflow 
rate. For example, it is believed that the energy radiated 
from shock waves provides a measure of the present driving 
power while the mechanical luminosity provides a record 
of the driving history. In this way, the ratio of emissions 
from vibrationally excited molecular hydrogen and rotation- 
ally excited carbon monoxide could deter mine the evolution- 
ary phase of t he protostellar exhausts iDavis et alJll993 
lYu et alJ koOO). We can determine here the validity of this 
statement, although these simulations are of a particular 
type and scale and so only begin to illustrate the full range 
of possibilities. For example, we take here a hydrodynamic 
highly-collimated jet and a uniform ambient medium that 
are both initially fully molecular. In the next section, we 
discuss more fully the characteristics of the jet simulations. 



2 ZEUS-3D WITH MOLECULAR COOLING 
AND CHEMISTRY 

2.1 The code 

We have modified the ZEUS-3D code, which usually updates 
variables with an explicit method, to include a semi-implicit 
method that calculates the molecular and atomic hydro- 
gen fr actions, according to the prescription of lSuttner et alJ 
(1997). We have also added a limited equilibrium C and O 
chemistry to calculate the CO, OH and H2O abundances 
jRosen fc Smithll2003h . Equilibrium CO and H 2 is a rea- 
sonable estimate at high density, consistent with the low nu- 
merical resolution of t he cooling layers behind shock waves 
jRosen fc Smithll2003D . 

The details of the many components of the cooling func- 
ti on and chemistry netwo rk are discussed in the appendices 
of lSmith fc Rosenl <|2003h . As an overview, we include cool- 
ing through rotational and vibrational transitions of H2, CO, 
and H2O, H2 dissociative cooling and reformation heating, 
gas-grain cooling/heating, and a time-independent atomic 
cooling function that includes non-equilibrium effects. We 
include an additional 10% by number of helium atoms, so 
the mean particle mass is 2.32 x 1(T 24 g. The dust temper- 
ature is fixed at 20 K. 

The ZE US code suffers from numerical errors. Recently, 
iFalld i2002T) has explored the appearance of rarefaction 
shocks which, as the name suggests, are particularly resis- 
tant and their influence cannot be removed by increasing 
the resolution. They may occur in locations where the local 
Courant number is near 0.5 and where the velocity is nega- 
tive. Inaccurate fluid patterns could result for adiabatic or 
MHD flows but isothermal flows are not significantly in er- 
ror. The present jet simulations are thus immune to these 
errors, with strong forward motion and strong cooling, espe- 
cially in the zones with relatively large Courant number (i.e. 
approaching 0.5). In any case, isothermal flows are an ap- 
proximation in which regions undergoing rarefaction are be- 



ing artificially heated (to counter expansion cooling). Cool- 
ing also generates structure through thermal instabilities on 
many scales which are not adequately featured in the numer- 
ical approximation. Nevertheless, we persist with the ZEUS 
code because of its high versatility, allowing the introduction 
of new physics, the discovery and removal of specific errors 
over many years, and its speed of execution. 

2.2 The fixed conditions 

Two models for the time evolution of the jet are studied. We 
designate these as Simulations D and I, with a decreasing 
and increasing jet density, respectively. 

In both 3D simulations, the computational grid spans 
2 x 10 17 cm along the jet axis (a;) and 2 x 10 1B cm in 
both transverse dimensions, with 1000 x 100 x 100 uniform 
zones in each direction. We initialise the jet with a nozzle 
radius, Kj, of 1.7 x 10 cm or 7.5 zones. The initial internal 
Mach number of the jets is 140 with a temperature for the 
fully molecular jets of 100 K; this combination gives a speed 
of 100 km s _1 . The ambient medium is molecular with a 
hydrogen nucleus density of 10 4 cm -3 . 

In addition, both jets are pulsed with an amplitude of 
± 30% and a period of 60 yr, with a sin(ui) dependence. A 
small precession with a 1° half- precession angle and a period 
of 50 years is included. A radial shear that reduces the speed 
to 70% of the core velocity at the jet edge is also included. 
The sinusoidal form of the pulsation and the parabolic f orm 
of the shear follow Eqns. 1 and 2 of IVolker et alJ lll999l) . 

2.3 The evolving parameters 

Two models for the time evolution of the jet are studied. We 
designate these as Simulations D and I, with a decreasing 
and increasing jet density, respectively. 

In Simulation D, we decrease the jet density abruptly, 
setting it to zero after a short period of time. We resorted to 
this dramatic variation of the flow because earlier attempts 
with more moderate decreases failed to show a significant 
deceleration. The end of the run may thus correspond to a 
highly evolved jet-driven outflow. Specifically, we allow the 
"jet" to propagate for a mere 20 years and then sever its 
ties to the life-giving momentum source at the inlet. The jet 
is initialised with a hydrogenic nucleon density of 10 s cm" 3 
(log rij — 5). In addition we have simulated other aborted 
jets, with flows terminated after 20 years but the flow was 
not initially precessed, and also a case where the flow was 
terminated after only 5 years and also not initially precessed. 

In Simulation I, we wish to show the consequences of 
increasing the mass density in a light jet flow that had al- 
ready propagated some distance. In this case, the density 
evolution follows a less abrupt change, determined by the 
following: 

Pi(t) = Pa +pitanh (— jj— ^) (!) 

where we set pi = 1.1484 x 10" 19 g cm -3 , p = 1.1716 x 
10~ 19 g cm -3 , ti — 800 years, and to = 50 years. This results 
in the logarithm of the jet hydrogenic number density rising 
from log rij — 3 to log rij = 5 over a 200 year period centered 
at a program time 800 years after the simulation has started. 
Given the axial grid length of 2 x 10 17 cm, this leads to the 
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Figure 1. Midplane cross sections of density. These are cross sections for z/Rj = 0.0 of Run D at t = 1300 yr (top) and Run I at t = 
1400 yr (bottom). Each is scaled to the same maximum and minimum and the scale, with darker shading indicating denser regions, is 
displayed below the cross sections. The complete computational domain in the two axes is displayed. In both panels, the vertical axis is 
the y-axis. 
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Figure 2. Axial cross sections of density from Run D (upper row) and from Run I (lower row). Each row is from the time displayed in 
the corresponding panel in Figure^ Each cross section is scaled to the same maximum and minimum and the scale, with darker shading 
to indicate denser regions, is displayed below the cross sections. We display the complete computational domain in the two axes shown. 
The distance between adjacent tick marks is 5 X 10 15 cm. The axial position in Rj is given in the upper left of each panel. The vertical 
axis is the j/-axis, and the horizontal axis is the z-axis, and so the view is from the jet inlet boundary looking down the jet axis (i.e., 
toward the +x-direction). 
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bow shock from the initially light jet being roughly half way 
across the grid when the mass flux is increased. 



3 INTERNAL STRUCTURE 

Slices of jet density down the midplane are presented in Fig- 
ure Both simulations feature a refocused leading region 
and wider trailing "shoulders" . We remark that the shoulder 
regions in both simulations shown in Fig.[Qhave progressed 
to only 7 x 10 16 cm, a small fraction of the jet leng th. In our 
study of non-evolving jets IIRosen fc Sm ith 2003), we were 
able to attribute the refocusing to the combination of both 
shear and strong molecular cooling, since simulations with 
atomic jets or molecular jets with no shear do not refocus. 
The ratio of the distance of the shoulders from the nozzle 
relative to that of the entire jet was found to be propor- 
tional to the jet-to-ambient density ratio. Hence, without a 
continuous strong power source in Simulation D, the bloated 
cocoon is stunted. 

Furthermore, the average speed of the leading bow 
shock, which almost traverses the grid in 1300-1400 years, 
is ~ 45 km s _1 . This is similar to a "light" jet simulatio n 
(jet-to-ambient density ratio of 0.1 in lRosen fc Smithl2 003). 
This average speed is still, however, nearly 90% faster than 
expected from ram pressure arguments due to the concen- 
tration of the jet thrust caused by the refocusing. 

From a plot of a time-averaged advance speed of the 
bow shock (Figure [3J, we see that the advance speed of the 
bow shock is not constant throughout either simulation. In 
Run D, the advance speed of the jet falls below 60 km s _1 
within the first 100 years, and follows a linear decrease in 
velocity with time until t ~ 800 years, where it continues 
a linear decrease but at a faster rate. The change of the 
deceleration rate is related to a broadening of the tip of 
the bow shock, which is at the leading edge of a region of 
refocused and accelerated flow after t ~ 200 years. At the 
end of Simulation D the bow shock is progressing at a speed 
of ~ 30 km s" 1 . This is an example of the ability of all of our 
aborted molecular jet simulations to maintain its momentum 
for the length of the computational grid (0.06 pc). Even in 
the case of the jet aborted after 5 years, we still found a 
speed of ~ 20-25 km s _1 after the 2000 year interval in 
which it crossed the grid. 

In contrast, the progress of the bow shock in Simulation 
I, as measured from the molecular emission maps presented 
below, increases early in the simulation. The speed initially 
drops to 20 km s _1 , but after only 50 years begins to rise. 
There is a significant increase at t ~ 300 years, after which 
the speed of advance of the bow shock is ~ 50 km s~ . The 
increase does not coincide with the onset of the higher jet 
density but, instead, coincides with the initial refocusing of 
the bow shock at its leading edge. Thus, the stripping of the 
low-speed jet sheath is an interesting means by which an 
acceleration in the proper motions of jet knots may result. 
Even after the jet density has increased, the average advance 
speed for the bow shock remains nearly constant until the 
simulation is terminated. 

The dense material in Simulation D is confined to the 
bow shock and swept-up shell while dense material is also 
associated with the internal shocks of the pulsed jet in Simu- 
lation I. The cocoon surrounding the jet behind the shoulder 




600 800 
t (yr) 



1000 1200 1400 



Figure 3. Time-averaged advance speed of the bow shock. The 
advance speed here is determined from the leading edge of the 
bow shock in integrations of H2 emission, similar to Figures HI 
and 151 



region has the lowest density, with the still-connected jet in 
Simulation I having a lower minimum density by a couple of 
orders of magnitude than the aborted jet in Run D. These 
overdense and underdense regions within the jet are also ap- 
parent in the axial cross sections of density in Figure|5| With 
the ability of the jet in Run I to generate a very low den- 
sity cocoon, the relatively small overdensity (still an order 
of magnitude) in the bow shock is difficult to see. Note twin 
concentric shell structures are formed in the cross section 
of Run I, visible in Fig. [5] at distances Kj = 54-91, where 
the powerful jet transmits invigorated shocks through the 



4 SIMULATED MOLECULAR LINE EMISSION 
4.1 The Spatial Distributions 

Images of line emission from molecules are shown for Simu- 
lation D in Fig. 0] and for Simulation I in Fig. 03 Integration 
is across the y-direction transverse to the jet axis, and the 
simulation time corresponds to the previous figures. The H2 
line emission is based o n a non-LTE approximation to the vi- 
brational populations ijHollenbach fc McKeelll979|). and the 
CO em ission is calculated from the formulae of lMcKee et alJ 
( 1982). In addition, we set a lower velocity limit of 1 km s" 1 
on the integration of CO emission lines in order to subtract 
out the undisturbed cloud. 

In FigureH we show the H 2 1-0 S(l) line, the CO R(l) 
and R(5) lines for the CO ground vibrational state (v = 0). 
In the H2 line, only the leading edge of the bow shock re- 
mains visible and, since Simulation D is relatively devoid of 
warm (T > 5000 K) gas, the image in the H 2 2-1 S(l) line 
is not shown. We display instead the CO R(5) line, which 
should highlight gas at intermediate temperatures (T ~ 100 
K) between the H2 line and the CO R(l) line (which empha- 
sizes the coldest, T ~ 5 K gas). The CO R(5) line image not 
only has an intermediate appearance but also shows some 
filamentary structure, which is brightest at the leading edge 
of both the bow shock and the shoulders. This relatively 




Figure 4. Appearance of Simulation D in three molecular emission lines. The luminosities from each zone are placed in bins in the 
viewing window that are the same size as the 3D zones used in the simulations (i.e., 2 X 10 14 cm). The axis being integrated is the 
2-axis, so the vertical axis in all the panels is the y-axis. 




Figure 5. Appearance of Simulation I in four molecular emission lines. See caption in Fig. [I] Since this simulation generated sufficient 
amounts of warm (T ~ 10,000 K) gas, we show the map of the H2 2-1 S(l) emission in the second panel. 
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warm gas is generated by the accelerating refocused region. 
In both of the simulations presented here, the moving CO 
R(l) shows the outline of the bow shock. 

Simulation I is sufficiently energetic to generate copious 
amounts of 10,000 K gas. Hence, in addition to the lines 
shown in Fig. 0] we also show the relative luminosity of 
the H2 2-1 line to the 1-0 line in Fig. |K| The emission in 
the H2 1-0 image is sparse and filamentary, showing several 
'minibows' especially at the edges of the shoulder region 
near x — 9 x 10 16 cm. Only a single pulse near the inlet 
is seen in the H2 1-0 map. Earlier pulses, farther down the 
jet, are easier to see in the ratio of the H2 lines, and are 
noticeable in the CO R(l) line, but very prominent in the 
map for the CO R(5) line. While the CO R(l) emission map 
does show a wider shoulder region well behind the leading 
edge of the bow, this region is much less pronounced than 
the appearance in the midplane density slice. 

4.2 Totals 

The evolution of the totals of the H2 emission lines and 
the CO R(l) line are shown in Fig. E] In Simulation D, a 
dramatic reduction in the H2 totals begins after t = 100 
yr. A flattening of the CO luminosity, which monotonically 
increases for non-evolving molecular jets, also occurs but not 
until t — 700 yr. Thus, a significant change in the ratio of 
the H2/CO total emission occurs shortly after the flow has 
been turned off. It is only when the H2 1-0/CO R(l) ratio 
approaches unity that the total emission from each emission 
line tends to flatten out. Before this occurs the ratio is a 
good age indicator for the aborted jet flow. 

The 80 year delay between the jet demise and the H2 
weakening is clearly the time to elapse before the leading 
bow stops being driven by the jet impact. The 700 year delay 
between the termination of the jet flow and the flattening 
of CO R(l) emission indicates that new ambient material 
continues to be snowploughed, as the momentum of a small 
fraction of high-speed gas is redistributed into a large frac- 
tion of low-speed gas. This is a true 'bullet' or ballistic phase. 
To test this, the deceleration of the bow shock in Run D as 
shown in Fig. [3] can be compared to the analytic motion of 
a bullet of speed Vf,, fixed mass M, cross-sectional area A, 
and initial density pto, moving into a medium of density p a , 
with a drag coefficient D c . The equation of motion is 



Alternatively, if the bullet is either disk shaped or does not 
expand at all, L is a constant and we find 



(2) 



dv b _ vl 
~dt ~ ~T 

where L — M/{D c p a A). Here L defines an effective column 
length of ambient gas required to slow down the bullet. First, 
we take ram pressure confinement of a bullet which main- 
tains its shape as the pressure decreases. This yields 
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where c 3 is the bullet sound speed and L is the initial value 
of L. Substitution and integration gives 
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For the particular conditions of Run D, L — (6 x 
10 16 /D c ) cm. To decelerate from 60 km s" 1 to 30 km s" 1 
in 1200 yr, then requires D c = 0.17-0.28, depending on the 
bullet model utilised. This drag coefficient is then consis- 
tent with the aerodynamic shape of the leading bow and 
thus shows that a large mass can be mildly disturbed by a 
high-speed jet. 

The duration of this phase is 35 times longer than the 
jet-driving phase, suggesting an effective jet column density 
which exceeds the deflected ambient column density by 35. 
We can roughly attribute this to the high jet density (10) 
and the refocusing (3.5). 

It should be remarked that decelerating bullets have 
been investigated in various astrophysical context s in- 
cluding Her big Haro objects {Norman fc Silkl Il979ft . ra- 
dio galaxies llSmith fe Normanlll980ft and planetary nebula 
(Pal mer et alJll996T) . Depending on how the bullet forms, 
stability is maintained through inertial confinement for a 
sound-crossing time or a shock-wave crossing time. There- 
fore, the high Mach number of the bullet in Run D permits 
a bullet to travel over 10 17 cm without great expansion. 

In Run I, periodic brightening in H2 luminosity occurs 
as the pulses in the flow pass through the shoulder region 
and energise the bow shock envelope. The increase in jet 
density that becomes significant for t > 750 yr shows up 
as an increase in the H2 output very soon thereafter (t = 
800 yr). No accompanying increase in the CO emission is 
noticeable. There is an unrelated early break in the rate of 
increase in the CO R(l) total luminosity, at t = 100 yr. The 
slightly increased rate is related to the shape of the envelope, 
which shows the effects of refocusing near the start of the 
simulation. 

The final integrated fluxes of the two simulations have 
contrasting behaviours. In Run D, the absolute luminosity 
in the 1-0 S(l) line of H2 is continued at the low level of 
L(l-0 S(l)) ~ 10" 5 L Q despite the lack of a jet. This com- 
pares to the value of ~ 10 -3 Lq for the eq uivalent simula- 
tion in which the jet power is maintained l|Rosen fc Smith! 
12003ft . In Run I, L(l-0 S(l)) reaches 3 ~ 10" 4 L Q and is set 
to increase. It is clear that a complex relationship exists be- 
tween the immediate jet power and L(l-0 S(l)). Whereas we 
found this ratio to lie in t he range 80-600 for non-evolving 
jets l|Rosen fc Smithl2003ft . we here find values between zero 
(Run D, after 20 yr) and 1.2 x 10 5 (Run I, 800 yr). 

The predicted values of L(l-0 S(l)) are compara- 
ble to those derive d in the unbiased statistical study of 
l|Stanke et alJl2002ft . where the majority of outflows possess 
L(l-0 S(l)) in the range 10" 4 -10" 3 Lq. Nevertheless, the 
absolute fluxes are quite low in comparison to a number of 
observations of well-known outflows. This is partly due to 
the compactness of the simulated outflows, limited by our 
need to follow molecule dissociation. We can argue that scal- 
ing up the spatial dimensions by a factor of 10 would scale 
the luminosities by a factor of 100. It remains to be shown, 
however, that the flow patterns would be maintained. 
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Figure 6. Integrated molecular emission in three lines for each simulation. The H2 1-0 emission is the solid line, the H2 2-1 is the 
dotted line, and the CO emission is the dashed line. 



Table 1. Mass Spectra Power-Law Dependences 



Sim. t(yr) 


type 


6 


range of log v 


7 


N a 


type 


e 


range of log v 


7 


N a 


D 1300 


mass 


15 


0.5-1.0 


1.19 


7 


CO 


15 


0.2-0.6/0.7-0.9 


1.27/1.92 


2/3 






30 


0.5-1.0 


1.17 


7 




30 


0.5-1.0/0.9-1.1 


1.17/2.74 


7/5 






15 


0.5-1.0 


1.19 


7 




45 


0.5-1.0/1.1-1.2 


1.22/3.58 


7/3 






60 


0.5-1.0 


1.10 


7 




60 


0.5-1.0/1.1-1.3 


1.16/3.07 


7/7 






90 


0.5-1.0 


1.13 


7 




90 


0.5-1.0/1.2-1.3 


1.21/3.18 


7/7 


I 1300 


mass 


15 


0.5-1.0 


1.05 


7 


CO 


15 


0.5-1.0 


1.20 


7 






30 


0.5-1.0 


1.01 


7 




30 


0.5-1.0 


1.22 


7 






15 


0.5-1.0 


0.99 


7 




45 


0.5-1.0 


1.22 


7 






60 


0.5-1.0 


1.01 


7 




60 


0.5-1.0 


1.21 


7 






90 


0.5-1.0 


1.03 


7 




90 


0.5-1.0 


1.21 


7 



N is the number of points in the velocity distribution, which was computed in 1 km s 1 bins, used to estimate the slope. 



5 VELOCITY DISTRIBUTION OF MASS AND 
MOLECULAR LINE EMISSION 

5.1 Mass "Spectra" 

CO velocity data for observed outflows are now invariably 
analysed by plotting the differential mass in velocity bins, 
with a few simple assumptions made to compute mass from 
CO intensity. One may classify how this mass decreases 
with velocity by fitting a power law, with the slope defined 
by dm/dv oc v K Two surveys containi ng many sources 
iYu et alJl2000l and lRidee fc MoordEoOll) conclude that a 
broken power law is often necessary, with a shallow slope of 
1.0-2.0 at low velocities and a steeper slope, in some cases 
7 is close to 10, at larger velocities. The velocity where the 
break occurs is usually ~ 10 km s _1 . One interpretation of 



this break speed is that it is the projected component of 
the critical shock speed, Vdcosd with Vd = 23 km s -2 (e.g., 
ISmithlll994h . above which all molecules are assumed to dis- 
sociate within the shock. 

We have performed a similar analysis for the velocity 
distributions of mass, of CO R(l) luminosity and of H2 1-0 
luminosity for each of the aborted jet simulations, including 
Simulation D, and for Simulation I. In the case of Simulation 
D, the analysis was performed at a program time of t = 
1300 yr, corresponding to other Simulation D illustrations 
presented here. In the increased-density jet of Simulation 
I, the tip of the jet has just barely reached the outer x- 
boundary at t = 1400 yr. We, therefore, chose to inspect 
data at t — 1300 yr. We list the results for a sample of 
viewing angles, defined as the angle of the jet axis out of the 
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plane of the sky and toward the observer, in Table We 
display the data for one viewing angle (15°) in Figure |7| 

Having only found low values of 7 in previous simu- 
lations involving molecular jets (although high values are 
possible for an atomic jet), we have chosen to investigate 
whether a reduced mass flux would give larger values for 
7. In Simulation D and the other aborted jet simulations, 
however, we still find shallow dependences of the mass dis- 
tribution on velocity i.e. small values for 7. The values for 
the aborted jet simulations are in general larger than their 
counterparts in the increasing jet-density simulation, with 
the values for both sets usually nearer 1.0 than 2.0. How- 
ever, large slopes are computed for high velocities in the only 
cases where a break could be seen, which are the CO-derived 
values for the aborted jet simulation. These results are not 
significantl y different t han those found for non-evolving jets 
jRosen fc Smithll2003ft . 

In the previous analysis for non-evolving molecular jets, 
the 7 determined for CO was nearly uniformly larger than 
that for mass. For the evolving jet simulations here, the de- 
pendence of CO luminosity on velocity is also usually steeper 
than the one for mass, but in some cases it is only equal 
(e.g., Run D at a viewing angle of 30°) or, for one aborted 
jet simulation, was lower. This occurred in a simulation of an 
aborted jet that propagated into an equal pressure ambient 
medium. The CO-derived 7 is lower than the mass-derived 
one at all viewing angles, and by as much as 0.6. However, 
in this simulation, the CO distribution can be better fit- 
ted with two power laws, with the higher velocity one much 
higher than the single linear fit to the mass distribution. 

Very little dependence of 7 on viewing angle is found 
for both of these simulatio n. Th i s is c ontrary to pre- 
vious jet simulations JSmith et a l\ |l997|. |Lee et alJ l200lL 
iDownes fc Cabrid I200A and iRosen fc Smith! l2003h . which 
have shown an inverse relationship between viewing angle 
and 7. However, the CO-derived slopes in the high veloc- 
ity range are again an exception. This group disagrees more 
strongly with the previous results and, for small viewing an- 
gles, 7 varies in direct proportion to the viewing angle. 

An excess at the highest velocities is found in Simu- 
lation I, as also found in previous studies of non-evolving 
molecular jets. This bump peaks at log v/km s _1 = 1.4 (for 
a viewing angle of 15°) and is associated with the internal 
shocks from the pulsed jet. In the aborted jet simulations, 
both the internal shocks and the high velocity bump are 
absent. 

Recent observations of H2 emission from pro- 
tostellar outflows have been similarly analysed 
JSalas fc Cruz-Gonzale3 120021) . As with the CO data, 
the H2 flux distributions show a break in behaviour, with 
the break between 2-17 km s _1 for the selection of sources. 
A flat or slightly rising flux is found at low velocities and a 
fall-off in flux with 7 in the range between 1.8 and 2.6 at 
high velocities. 

From Fig. [7] we see that Run I may best match these 
observational constraints. However, the H2 flux distribution 
in Run I rises one order of magnitude in brightness from the 
first H2 datum in the figure to where it turns over near v 
= 10 km s -1 . This is much larger than that found in the 
observations, although it is a smaller rise than found in the 
non-evolving molecular jets simulations. Additionally, the 
best fit line between logw = 1.0 and 1.3 has 7 = -1.6, which 
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Figure 10. Velocity channel map for Simulation D in CO R(l) 
emission. The jet is viewed 15° out of the plane of the sky toward 
the observer. The central radial velocity in km s — 1 is indicated on 
the left side of each panel. Four orders of magnitude in luminosity 
are displayed in each panel. The log of the maximum luminosity 
in ergs s — 1 is shown on the right side of each panel. Each bin 
spans 2 X 10 14 cm in x' , the jet axis projected onto the plane of 
the sky, and in y. 



is in line with the observations. This could be an adequate 
explanation for sources viewed with their jet axis near the 
plane of the sky, such as HH 212. However, at larger viewing 
angles (Fig. |7| shows a viewing angle of 15° only), the low- 
velocity rise steepens. In addition, at intermediate velocities 
the comparison becomes less good as the viewing angle is 
increased. C-type shock physics may solve these problems 
since ambipolar diffusion would generate considerably more 
emission near zero radial velocity in the absence of discon- 
tinuous acceleration. The effect this has on the global prop- 
erties of an outflow remain to be calculated. 



5.2 Position- Velocity Maps 

We display position-velocity maps in H2 and CO emission 
lines for Simulation D (Fig. |HJ and Simulation I (Fig. 0. 
The H2 1-0 position-velocity map for Simulation D shows 
emission only at the leading edge of the jet. This emission is 
spread over 10 km s" 1 and 1 x 10 16 cm along the projected 
jet axis and is characterised by two peak s. Note that steady 
state bow shock models (see Fig. 14b of ISmith et al.ll2003l) 
predict two peaks, which correspond to (1) the projected 
edge-on leading edge near zero radial velocity and (2) the 
blue-shifted emission from the warm material projected well 
behind the bow apex. 

The pulsed jet in Simulation I shows intermittent H2 

1- and 2-1 emission along its length, with many emission 
sites spreading out over a large range of velocities. The H2 

2- 1 map shows that the sites of bright H2 1-0 emission are 
also bright in the more highly excited line. 

The emission in the CO lines shows a pair of "Hubble- 
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Figure 7. Distribution of velocities into bins of mass and two molecular line luminosities for a viewing angle of 15°. Displayed are 
the distribution of mass (top panel), CO luminosity (middle panel) and H2 1-0 emission (lower panel). Each velocity bin is 1 km s _1 
wide. The designation for the data presented in each panel is shown in the top panel. Each run may be represented twice within each 
panel, with data from both positive and negative radial velocities included. Naturally, the smaller ranged data is for the positive radial 
velocities (which could contribute to a "red" lobe, while the other is from a "blue" lobe). 
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Figure 9. Position velocity map for Run I in four molecular emission lines. The jet axis is 15° out of plane of the sky toward the 
observer. Each bin spans 2 X 10 14 cm in x', the jet axis projected on the plane of the sky, and 1 km s — 1 in radial velocity. 



law" like dependences for Simulation D, but disconnected 
from the nozzle. The outer envelope of the first increases 
relatively quickly with position, and is followed by a shal- 
lower increase for much of the remainder of the jet. While the 
emission in the CO R(l) line is nearly uniform for non-zero 
velocities, there is some structure in the CO R(5) position- 
velocity diagram. The bright spot near x' = 7 x 10 16 cm 
is from the advancing edge of the shoulder region seen in 
Fig. 2] and the peak at x = 5 x 10 16 cm is from the trail- 
ing filamentary structure. Both of these bright features have 
radial speeds near zero. 

For Simulation I, the CO R(5) diagram shows more 
structure than the diagram for CO R(l). Both show a promi- 
nent series of triangular shapes that distinguish the Hubble 
Law behaviour of a pulsed jet. There is evidence for a con- 
nective loop of emission between successive Hubble Law re- 
gions, until the triangles merge downstream, filling in all ve- 
locities, and the Hubble Law is then hidden. The maximum 
luminosity in the CO R(5) map falls in between those of the 
H2 1-0 and the CO R(l) maps. Naturally, all of the maxima 
are higher than the corresponding values in the aborted jet 
simulation. 



5.3 Velocity Channel Maps 

We display CO R(l) velocity channel maps for Simulation D 
in Fig. 1101 and Simulation I in Fig. Illl with a resolution of 1 
km s~ x . While the CO emission from Simulation D does not 
span a large velocity range, the morphology does change as 
one moves from relatively high radial velocity (10 km s _1 ) to 
low radial velocity. At high radial velocities, only the leading 
edge of material is visible, and as the velocity is decreased 
more emission appears closer to the inlet. This is another as- 
pect of the Hubble Law profile seen in the position-velocity 
map. Near zero velocity the entire envelope of the jet and 
cocoon of shock ambient material is represented, along with 



the shoulder region seen in the density cross sections and in- 
tegrated CO maps. This defines a T-type morphology. Even 
the filamentary structure seen in the integrated emission 
map for CO R(5) is seen to some extent in these CO R(l) 
velocity channel maps. 

The CO R(l) velocity channel maps for Simulation I re- 
veal morphologies in the different vel ocity ranges more typ i- 
cal of our previous non-evolving jets llRosen fc Smithl2003fl . 
At high radial velocities, only the internal shocks from the 
pulsed jet are seen and at low velocities the entire envelope 
of jet plus cocoon material is present. At low velocities, the 
gap between the nozzle and the onset of emission decreases 
and a V-type global morphology is identifiable. At interme- 
diate velocities, a spinal structure is present. 

The twin protuberances at the front of the jet are from 
the differing advance speeds associated with the extrema of 
the evolving jet density in Simulation I. As confirmed from a 
sequence of integrated H2 1-0 emission line maps (Fig. 1121 . 
a faster-moving bow shock leaves the path of the earlier 
slower-moving initial bow shock. This departure occurs at a 
fairly late time in the evolution of the jet (t = 1300 yr, or 550 
yr after the onset of the increase in jet density at the inlet) 
and well into the refocused leading part of the bow (x = 1.6 
x 10 17 cm). The average speed of the new bow is therefore 
close to 100 km s" 1 . An average speed for the propagation 
of the new bow shock equal to that of the (time-averaged) 
jet speed is not unreasonable in this simulation with the 
very light initial jet initially evacuating the path and with 
an acceleration from the refocused region. 



APPLICATION: 
DIAGNOSTIC? 



AN EVOLUTIONARY 



The structures and features detailed above can provide clues 
to distinguish the evolutionary phase of a specific outflow. 
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Figure 11. Velocity channel map for Simulation I in CO R(l) 
emission. The jet is viewed 15° out of the plane of the sky toward 
the observer. The central radial velocity in km s _1 is indicated on 
the left side of each panel. Four orders of magnitude in luminosity 
are displayed in each panel. The log of the maximum luminosity 
in ergs s _1 is shown on the right side of each panel. Each bin 
spans 2 X 10 14 cm in x' , the jet axis projected onto the plane of 
the sky, and in y. 
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Figure 12. Birfucating bow shocks shown in the H 2 1—0 line for 
Simulation I. We display the portion of the computational grid 
with x = 1.4-2.0 X 10 lr cm. Adjacent tick marks are 1 X 10 16 
cm apart. An arrow points out the progress of the bow shock 
associated with the denser flow. 



Statistically, however, the integrated properties can pro- 
vide the most employable measures. In particular, the ra- 
tio L(H 2 1-0 S(l))/L(CO R(l)) can be determined from 
simulations and compared to the observationally derived 
L(H2)/L mec h. The total molecular hydro gen emission L(H 2 ) 
can be estimated from L(H 2 1-0 S(l)) JSmithlll993) . pro- 
vided the near-infrared extinction can be determined. The 
mechanical luminosity for both lobes in a bipolar outflow is 
derived from formulae of the form L mec h = M co v^ /(2l) 
where M co is derived from L(CO R(l)), assuming a CO 
abundance and temperature, v co represents the appropri- 
ately averaged speed of the mass that dominates the outflow 
energy (and must account for the source orientation when 
derived) and I is the outflow size. 

From the ratio of mass to L co in the log v bins of Fig- 
ure [7| we estimate from the simulations that Mco for one 
side of a bipolar outflow ~ 300 L(CO R(l)) in solar units. 
Furthermore, given the consistently low values of 7, we take 
v co = 20 km s" 1 . Finally, we fix L(H 2 1-0 S(1)/L(H 2 ) = 
0.05. It must be remarked that a rather high number of 
assumptions precede a comparison of simulation and ob- 
servation. Nevertheless, it is clear that the ratio L(H 2 1— 
S(l))/L(CO R(l)) ~ 3-8 for Run D after 1000 yr, while 
L(H 2 1-0 S(l))/L(CO R(l)) ~ 30-40 for Run I after 1000 
yr. Given the above factors these convert to: L(H 2 )/L mec h 
~ 0.02-0.05 for Run D and ~ 0.2-0.3 for Run I. 

We confirm that strong H 2 flows are thus associated 
with Run I. For example, we find L(H 2 )/L mec ; l ~ 0.3 
for Cepheus E, on e of the dynamically youngest outflows 
jSmith et all2003l) . Moreover, the H 2 emission is widespread 
and the CO position-velocity diagram displays structure cut- 
ting through all radial velocities at all positions, rather tha n 
the triangular Hubble-law features llLadd fc Hodaprj|ll997f) . 
consistent with the Run I predictions. Finally, we note that 
the CO 7 values for Cepheus E are 0.99 and 1.8, for the 
two lobes, although these va lues are quite sen sitive to the 
chosen local standard of rest jSmith et alJll99H) . Hence, ob- 
servations of Cepheus E are consistent with our results for 
a newly-forming outflow. 
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To date, CO emission has not been fully re- 
solved in many outflows. Velocity channel maps are thus 
rare. A highly collimated CO outflow emanates from 
IRAS 21391+5802, with a V-shaped m orphology on veloc- 
ity channel images feeltran et alJ l2002l. The V-shape sug- 
gests that it is also an immature outflow in which the central 
source is increasing in power. We would thus expect to find a 
strong H2 outflow, as is, indeed, the case: the total shocked 
emission is estimated to be 4Lq whereas the mechanical 
powe r estimates lie in the range 0.15-1.2 L© iNisini et alJ 

In contrast, there are many examples of outflows 
with very weak total H2 emission. These include the 
HH366/B5IRS1 outflow in the Perseus cloud complex, 
which possesses high collimation. Hubble-law behaviour, 
consistent with a single outburst in the past is also present 
<Yu et alJll999T) . The H outflow in OMC-2 is another exam- 
ple, where the presently-observed H2 jet does not appear to 
be the driver for the CO outflow <Yu et alJl200oh . 

IRAS 18148-0440 in the L483 core is thought to be a 
protostar undergoing the transition from Class to Class 
1. Its CO outflow may possess a T-shaped morphology, al- 
though background cloud subtraction makes this uncertain 
jTafalla et alj EoOO). It belongs to a class of objects charac- 
terised by low velocity CO and no C O bullets, as also p re- 
dicted here for a transitional outflow jTafalla et alJl2000l) . 



7 SUMMARY 

We have reported on two jet simulations that model drastic 
jet density evolution over tens or hundreds of years. The ba- 
sic physical structure is similar to th at previously found fo r 
non-evolving sheared molecular jets llRosen fc Smitbl boOS). 
with an accelerated leading bow shock followed by a wider, 
slower shoulder-like region. More specifically, both of the 
simulations on display here resemble a light molecular jet 
simulation. The primary reason for this is that the shoulder 
region makes only slow progress relative to the leading edge 
of the bow shock. 

We have explored predicted features relevant to high 
resolution near-infrared (H2 rovibrational at 2.12/im emis- 
sion), far-infrared (CO R(5) rotational emission at 435/^m) 
and millimetre (CO R(l) rotational emission at 1.3mm) 
wavelength observations. As expected, the H2 lines reveal 
patchy details of the structure and the CO lines show the 
general outline of material encompassed by the bow. CO 
line emission from mid-level rotational states show both the 
general outline and the internal shock structure. Owing to 
the small amount of higher temperature gas in the aborted 
jet simulation, the H2 emission is limited to the very tip of 
the accelerated bow shock. The predicted values of the 1- 
S(l) 2.12/im luminosity are comparable to those derived 
fo r the majority of out flows in the unbiased statistical study 
of IStanke et"alN2002l) . 

The evolution in the integrated luminosity in these 
molecular lines is evident in both simulations. The changes 
in the integrated values of H2 occur at program times that 
are within 100 yr after the modification of the inflow jet 
density. Of particular note is the ratio of the integrated H2 
to the integrated CO R(l) line in the aborted jet case, which 
decreases dramatically with time. Thus, in principle this ra- 



tio can be used to estimate the age of an outburst or aborted 
molecular flow from a protostellar core. 

We have also computed velocity distributions in mass, 
CO luminosity, and H2 luminosity. Here the main results 
are: 

1. The slopes of these distributions are small, closer to 
one than two in most cases. Steeper slopes (7 = 2-3) are 
seen at high velocities in CO in the aborted jet simulation. 

2. The slopes for the aborted jet case are larger than 
those for the increasing jet density case, although not by a 
large amount for the simulations list ed in T able U 

3. Similar to non-evolving jets dRosen fc Smitbll2003t) . 
7co is usually greater than 7 maS s, but there are significant 
exceptions. 

4. In the aborted jet simulation, there is no excess in 
the amount of mass or CO at the highest velocities; we have 
seen this excess that is associated with the internal shocks 
of a pulsed jet. 

The position-velocity maps for the aborted jet case are 
relatively featureless. In CO emission, a triangular Hubble- 
Law outline is filled in by the emission but does not include 
the low velocity emission near the jet inflow boundary. The 
CO R(5) position-velocity map shows some fine structure, 
with brighter regions associated with the warmer parts of the 
flow: the leading part of the bow shock and some internal 
structure near the front of the shoulder region. 

In the increasing jet density simulation, the position- 
velocity CO maps show the series of Hubble Law regions 
associated with successive pulses in the jet. Additionally, 
the CO R(5) map has more structure than the CO R(l) 
map. In the H2 position-velocity map, there are numerous 
features that are narrow along the jet axis but cover a large 
range of velocities (the vertical lines in Fig. |5J. 

The CO R(l) velocity channel maps for the aborted jet 
simulation have emission confined to a small range of veloci- 
ties. The sequence of channel images shows a T-morphology 
at low velocities and the characteristics typical of a Hubble 
Law at high velocities. The channel maps for the increasing 
jet density simulation possess a V-morphology, with the in- 
clusion of an additional protuberance from the new denser 
flow that carves a fresh path near the previous flow. 

These features, combined with the levels of inte- 
grated emissions, suggest that jets in certain outflows (e.g. 
IRAS 21391+5802 and Cepheus E, see Sect.[5J are increasing 
in power, while other jets (e.g. that from IRAS 18148-0440 
in L483) are on the decline. 

Not many detailed velocity channel maps have as yet 
reached the literature. This will change shortly, and with 
the advent of new telescopes such as SIRTF, FIRST and 
ALMA, imaging and spectroscopy will reveal the properties 
on fine scales in many molecular lines of intermediate and 
low excitation, such as in the CO R(5) and R(l) lines. Here 
we have determined how these properties can help determine 
the jet evolutionary stage. 
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